Simulations of Spinodal Nucleation in Systems with Elastic Interactions 
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Systems with long-range interactions quenched into a metastable state near the pseudospinodal 
exhibit nucleation that is qualitatively different than the classical nucleation observed near the 
coexistence curve. We have observed nucleation droplets in our Langevin simulations of a two- 
dimensional model of martensitic transformations and have determined that the structure of the 
nucleating droplet differs from the stable martensite structure. Our results, together with experi- 
mental measurements of the phonon dispersion curve, allow us to predict the nature of the droplet. 
These results have implications for nucleation in many solid-solid transitions and the structure of 
the final state. 
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An outstanding challenge in the physics of solid-solid 
phase transformations is to understand the nucleation 
and growth of the transformed phase in strain based ma- 
terials such as martensites, ferroelectrics, and multifer- 
roics. Although ideas based on classical nucleation the- 
ory have been invoked to describe nucleation phenom- 
ena in these materials 

B S H II II 0, it is possible 
that the presence of long-range elastic forces substan- 
tially influences the probability of nucleation and sub- 
sequent growth that determines processing and material 
behavior. For example, the transformation susceptibil- 
ity determined by nucleation is one of the basic factors 
for the hardenability of steels. Even though heteroge- 
neous nucleation, which is sensitive to the distribution 
of appropriate structures in the parent phase from which 
product phase nuclei may be triggered, has long been 
considered important for martensites 0], homogeneous 
nucleation also has been observed if the transformational 
driving force is sufficiently large. However, due to the dif- 
ficulty of experimentally determining nucleation droplet 
structure in martensites Qf, very little is quantitatively 
known about the morphology of the nucleating droplet 
or "embryo," including its size distribution or nucleation 
rate for homogeneous or heterogeneous nucleation. The 
aim of this work is to probe for the first time the nature 
of the nucleating droplet and associated fluctuations by 
performing mesoscale simulations using realistic nonlin- 
ear models for martensites. We find that the classical 
theory does not accurately describe the structure of the 
nucleating droplet, but that concepts associated with nu- 
cleation near a spinodal 0, account for the observed 



droplet morphology. In addition, our work provides the 
basis for investigating nucleation in systems where long- 
range elastic forces crucially determine the morphology. 



Long-range elastic interactions are important in strain 
based materials, such as martensites, and result from the 
requirement of compatibility of strain components which 
is necessary to preserve the continuity of the elastic me- 
dia [Tol ITU . We refer to systems with long-range interac- 
tions as near-mean-field, and use several characteristics 
of mean-field theory to study nucleation in these sys- 
tems. Mean-field systems have a well-defined spinodal, 
the limit of metastability; near-mean-field systems have 
a pseudospinodal, which becomes better defined as the 
range of interaction increases • 



Spinodal nucleation, that is, nucleation close to the 
pseudospinodal, is predicted to produce a ramified 
droplet with a small amplitude H, Q . The droplet need 
not have the same structure as the stable phase 0, , 
unlike in classical nucleation where the droplet is com- 
pact and has the same structure as the stable phase ■ 
Spinodal nucleation has been observed in molecular dy- 
namic simulations of simple models [T^ llT j , but the con- 
cepts have never been tested on more realistic represen- 
tations of materials such as martensites. 
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We model a martensitic transformation using a 

Ginzburg-Landau (GL) free energy of the form [lOl llll 
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where is the deviatoric strain, ei is the shear strain, 62 
is the compressional strain, U{p = r — f*) is the Fourier 
transform of 



U{k) 
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T = {6 — 6c) /{9q — 9c), 9 is the dimensionless temper- 
ature, 9c is the critical temperature where the = 
austenite minimum of Fq disappears, 9o is the tempera- 
ture where the three minima of Fg are degenerate, and R 
is the range of the interaction. For metals R is quite large. 
The quantity Fes can be written as nonlocal surface and 
bulk terms in (j){r) by using the St. Venant compatibility 
equations ^3 ■ If the width of the interface of the droplet 
scales as the correlation len gth f , the surface term can 
be neglected for ^ >> i? [T3lll8l |. The exponential cutoff 
has been added to U{p) in Fes to simulate defects. All of 
the variables are dimensionless and scaled, as described 
in detail below. 

We use overdamped Langevin dynamics so the equa- 
tion of motion for 0(r, t) is 



d^{r,i) ^ dF[^{r,t)] 
dt S(t>{f, t) 



+ ar,t), 



(3) 



where the Gaussian noise ( is related to the dimensionless 
temperature 9 by the fluctuation-dissipation relation. 
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Several predictions have been made [I^ about the crit- 
ical droplet when the system is quenched to just above 
the spinodal temperature, Tg, so that At = t — Ts « 1. 
In d = 2 the nucleation barrier is proportional to the 
Ginzburg parameter G ~ C^0^/(^'x) ~ i'^ 4>'^ / {9 cX) » 1, 
where x ~ 2/ At is the susceptibility. Because w 
•y/ At/4 << 1, the droplet is difficult to distinguish from 
the metastable background. For a > 0, ^ « y^a/ (4At) 
a nd T, = Ta^, so G « aT/{89c). If a < 0, ^ w 
^|a|/(2AT) and t, = |a|V(86), so G « |a| At/(46'c). 
The droplet is predicted to be modulated at the largest 
(real) value of the wavenumber fcp at which the structure 
function S{k) w [t ± \a\P/4: + bk^/8 + Uik)]-^ diverges. 



If a > 0, fco = and the droplet is homogeneous, with 
(f>{r) in the droplet either everywhere positive or nega- 
tive. If a < 0, fco = vH/^) ^i^d the droplet is modulated 
with wavelength w = 2n^ybJ\a\ << ^, with alternating 
regions of positive and negative (t>{r). 

One of our main goals is to simulate martensites as 
realistically as possible. To this end, we need to relate 
the dimensionless simulation parameters in Eq. ^ to 
empirically accessible parameters. By extending the po- 
tential uT] to include a second gradient term, the three- 
dimensional (elastic) free energy near the critical temper- 
ature in terms of measurable quantities is 
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where oq is the crystal lattice spacing, L is the linear di- 
mension of the system, and Gi through Gg, ki/oq, and 
K2/aQ are elastic constants in units of N/m^. L/oq >> ^ 
to minimize finite size effects. All the parameters in 
Eq. ((SJ can be determined empirically [Tlj : for example, 
Ki/aQ is determined from the curvature at small k of the 
phonon dispersion curve obtained from neutron scatter- 
ing experiments . We scale all the elastic constants by 
Aq = 9G|/(128G|) and the strains by eo = y/^C^J{AC^ 
to make the two martensite minima of the homogeneous 
part of F in Eq. (|lb|l near unity, and we scale all distances 
by flQ. We define the dimensionless variables r = r/ao, 



L/ao, F = F/iAoal), ei 



C3eg/(2Ao) 



ei,2/eo, = ea/eo, 
a = 2Kiel/{Aoal), 



L 

Al,2 = 0x^2^1/ Aq, t = 

and b = 2K2eo/(^oao)- 

To find the proper scaling for the time, noise, and the 
temperature, we use the Langevin equation for e3(r, i), 
which is Eq. ||2Jl with £3, F[(j)] -^[£3], and 
an explicit friction coefficient 7 multiplying the time 
derivative. The dimensionless time and noise are t = 
Aoaot/{LjeQ) and ( = Leo(/{Aoao). The factor of L 
is a result of going from ez{x,y,z,t) to (j){x,y,t). The 
fluctuation-dissipation relation Eq. Q requires that we 
define 9 — kBT/{AoaQ). 9c and are defined similarly 
with T replaced by Tc, the critical temperature, and Tq is 
the temperature at which the three minima of the homo- 
geneous part of F in Eq. jSJl are degenerate. Note that 
we can change the effective temperature of the simula- 
tions either by changing T or ag. We drop all tildes in 
the following. 

To discretize the Langevin equation, we used a sim- 
ple forward Euler method j23| for the time derivative. 
Higher order algorithms take more time and give similar 
results. The random noise ^ is computed by multiply- 
ing the standard deviation of the noise by a random 
number chosen from a Gaussian distribution with unit 
variance, Gij^a On the lattice 5^{f— r') « l/fa^ 

and S(t — t') « ^/St, so C(^i^) ~ <^cGij,a, where cr^ = 
\J29 1 (Sx'^St). The treatment of the spatial derivatives 
and the nonlocal term is not as simple, however. 
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The nonlocal term in the Langevin equation that arises 
from Eq. (le) is a convolution with the kernel K{k) — 
J(P pU{p)e~^P^^^e'^^'''P. Here K{k) needs to be com- 
puted only once. The exponential factor e^'^'/^ is com- 
puted with the origin chosen at the center of the lattice 
so that it obeys periodic boundary conditions. On the 
discrete lattice the indices i and j correspond to x and y, 
I and m correspond to kx and ky, and a corresponds to t. 
We write (j){x,y,t) ^ U{kx, ky) Ui^^, e'l'^l/^ 

Eij, and K{kx,ky) ki^^n ~ T[Ei^-iT^^{Ui^,n)], where 
J- and J^~^ represent fast Fourier transforms. The con- 
volution integral is computed as: 

/ dV'{/(p)e-|p1/«'0(r-',i) « (6) 
Jo 

The accuracy of our algorithm depends on treating 
both spatial derivatives to fourth order in 6x, taking 
care to include cross terms in V'*^. If we define (s'fj — 

(2) 

4>i+l,3 + (t>i-l.,j + (f'i.j+l + (t'l.j-l, ©ij = 4'i+2,j + + 

(t'i-2,j+4>i,j-2, and (8)ij = (j)i+i^j+i+(j}i+ij-i+4>t-i^j+i + 
then 

" 12^^- ®S +16 -600,,] (7a) 
V^0.,, « J^i®S - 8 ®^^} +2 +20(/..,,]. (7b) 

We use the values = 268 K, Tq = 290 K, and 
A2 = 2A-\ for all our simulations; these values correspond 
toFePd 111 so that our simulations are as realistic as pos- 
sible. Wc also used L — 64 and 6x — 0.5. Our procedure 
is to fix Sx and St = 0.01 and choose the values of the gra- 
dient parameters a and b so that the numerical solution 
is stable. Our numerical solutions were checked for ac- 
curacy by comparing the simulation results to the exact 
analytical solution for the linear case without the noise 
and nonlocal terms. Numerical stability for the complete 
equation of motion was checked by varying St and Sx. 
We limit ourselves to a and b such that the Ginzburg 
parameter G ~ 5 and ^ w 16 = L/4. We choose Ai < 1 
so that the core of the droplet will be more visible. We 
varied the nucleation rate primarily by changing ap. 

For a > we take a = 6.32, b = 0.01, R = 6.4, r = 
6.17 X 10-3, Ai = lx 10-3, and ao = 2.1544 x lO-^m. 
FePd has a crystal lattice spacing of « 3 A, so this value 
of ao corresponds to a coarse-graining factor of 70 and 
sample size of about 1.38 ^m. The value of t corresponds 
to T = 268.14 K, and a corresponds to ki/oq = 2.79 x 
lO^N/ni . The value of a is reasonable compared to the 
value of 2.5 x lO^^N/m^ for FePd quoted in Ref. 'Tl'|. 

For a < we take a = -1.61, b = 0.652, i? = 4, 
T = 0.5, Ai = 0.6, and ao = 2.7 x 10-^ m. These values 
correspond to the modulation wavelength it; w 4, = 
0.49695, T = 279 K, and At = 3.048 x 10-^. 

We begin the simulation with 0(r) chosen at random 
around the metastable austenite minimum at = 0. 



Within a short time, (j){f) equilibrates at the chosen tem- 
perature. On a time scale much longer than the initial 
equilibration, the thermal noise causes a critical droplet 
to appear. By looking at the evolution of the spatial 
average, <(jp(t)>, we can make an estimate of the nu- 
cleation time. From Fig. ^ we see that the system went 
from a metastable state with <(jp> ^ to the stable 
state with <0^> 1 at the time t = 2430. Because we 
expect the amplitude of the droplet to be close to that of 
the metastable phase, we expect the nucleation time to 
be before the rapid increase in <(j>'^>. 




2360 2380 2400 2420 2440 2460 2480 2500 2520 2540 

Time 



FIG. 1: Evolution of <(f>'^{t)> for a > 0; the original run and 
three interventions at t — 2355. <4>'^{t)> of the original run 
grows rapidly from its metastable value to the stable value at 
t ~ 2430. If the intervention run is similar to the original run, 
the same droplet is assumed to have grown; otherwise, it is 
assumed that the same droplet did not grow. Some interven- 
tions are ambiguous. In this run the intervention time of 2355 
is close to the estimated nucleation time of — 2356 ± 35. 

To determine the nucleation time more precisely, we 
use an intervention technique |2l|. Because the nucle- 
ation droplet is a saddle point |^, the droplet has 
an equal probability of growing to the stable state or 
shrinking back to the metastable state if we perturb the 
system at the nucleation time. We can use this saddle 
point property to help us find the nucleation time. We 
restart the simulation at the estimated nucleation time, 
and integrate the equations of motion using a new se- 
quence of random numbers for the thermal noise. Our 
criterion is that if 8 ± 4 of the 16 runs with different ran- 
dom number sequences show that <(jP> grows at roughly 
the same time as in the original run, then the interven- 
tion time is equal to the nucleation time. We can then 
look at snapshots of (f){x, y) to see if we can identify the 
droplet. 

For a > 0, the nucleation time is t+ = 2356 ± 35, 
corresponding to Fig.^ A snapshot of 0(a;, y) &it = 2356 
is shown in Fig. 2a. Even though only the core of the 
droplet is visible above the noise, we can see that the 
droplet is homogeneous, as predicted. 

For a < 0, the nucleation time is estimated to be t~ — 
1236.5 ±1.5. The error in the nucleation time is much 
less than for a > because the slower growth in <(/)^> 
for a < makes fewer of the interventions ambiguous. A 
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FIG. 2: (Color) From left; (a) The deviatoric strain 4>{x,y) for a > at the estimated nucleation time — 2356. The 
amplitude of ranges from —0.08 (dark blue) to 0.025 (yellow). The part of the droplet visible above the noise is clearly 
homogeneous, showing no modulations between the two low temperature minima (red and blue), (b) The deviatoric strain 
4>{x, y) for a < at the estimated nucleation time t~ = 1236.5. The amplitude of (j> ranges from —0.2 (dark blue) to 0.2 
(red). In the part of the droplet that is visible, the modulations between the low temperature minima (red and blue) have the 
wavelength predicted by Ref. [T^ to be w ~ 4. 



snapshot of the field &tt — 1236.5 is shown in Fig. 2b. As 
predicted in Ref. the droplet has modulations with 
wavelength w « 4. This modulation is different from the 
twinning in the stable phase which occurs at wavelength 
A - ^/I= 8 [13. 

These results are important for nucleation in many 
functional materials with long-range interaction in which 
strain couples to some other physical variable, e.g. polar- 
ization (in ferroelectrics), magnetization (in magnetoe- 
lastics), and other multiferroics. We found that in elastic 
systems the description of nucleation is subtle due to the 
presence of bulk/interface elastic compatibility contraints 
that are manifested as long-range interactions. In sum- 
mary, when the austenite is quenched to near the pseu- 
dospinodal, the structure of the nucleation droplet is dif- 
ferent than the structure of the stable martensite phase. 
If the curvature of the phonon dispersion curve at small 
k is positive, then a > and the droplet is homogeneous. 



If the curvature is negative, then a < and the droplet 
is modulated with a wavelength w. The parameters we 
used are consistent with austenite to martensite transi- 
tions in FePd. We conclude, therefore that the classical 
nucleation picture is not applicable to these transitions 
and that the spinodal nucleation scenario is the better 
approach to understanding these transitions. 
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